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Introduction 


The  main  objective  of  this  project  is  the  development  of  accurate  physical 
models  and  efficient  numerical  algorithms  suitable  for  diffraction  tomographic 
reconstruction  of  the  compressibility,  acoustic  attenuation,  and  mass  density 
of  the  prostate  through  a  multielement  transurethral  ultrasound  transceiver. 
Unlike  conventional  ultrasound  imaging,  which  is  non-quantitative  and  af¬ 
fected  by  speckle  artifacts,  three-dimensional  tomographic  solution  of  the 
inverse  acoustic  scattering  problem  has  the  potential  to  quantitatively  recon¬ 
struct  the  detailed  acoustic  properties  of  the  prostate  from  measurements  of 
the  scattered  radiation  field.  However,  numerous  challenges  must  be  con¬ 
fronted  before  such  a  proposition  becomes  feasible  in  real-world  applications. 
It  is  well-known  that  problems  in  inverse  scattering,  in  addition  to  being 
mathematically  complex  and  computationally  intensive,  are  also  particularly 
ill-posed  and  ill-conditioned. [1,  2] 

In  the  weak  scattering  limit,  such  problems  are  tractable  within  the  so- 
called  Born  approximation,  in  which  the  pressure  field  internal  to  the  scat- 
terer  is  approximated  by  the  pressure  field  of  the  incident  wave.  This  approxi¬ 
mation  assumes  negligible  loss  of  the  incident  field  energy  through  scattering, 
which  is  effectively  the  single  scattering  limit  where  7  — >  0,  and  \/L  »  1. 
The  Born  approximation  is  clearly  inappropriate  for  imaging  an  object  such 
as  the  prostate  which  may  not  be  a  weak  scatterer  at  all,  and  for  which 
multiple  scattering  must  be  taken  into  account.  For  this  reason,  we  have 
taken  an  approach  involving  solution  of  the  full  multiple-scattering  form  of 
the  nonuniform  wave  equation; 

7K(r,  ()  +  V  .  (7,(r,  «)Vp)  -  2V  •  M  •  V,  (1) 
recast  in  the  form  of  a  Lippmann-Schwinger  integral  equation: 

Pa;(r)  =  Pmc(r)  +  /  (A:^7KPa;G'a>(r|ro)  +  7pVoPa; '  VoG<^(r|ro))  dro,  (2) 
Jn 

as  discussed  in  Morse  and  Ingard.[3]  Our  initial  implementation  has  only 
considered  variations  in  compressibility  in  the  scatterer,  neglecting  the  second 
term  in  Eq.  2;  incorporation  of  density  inhomogeneities  in  the  forward  model 
will  constitute  a  significant  component  of  work  in  the  upcoming  year. 

The  characteristic  physical  scales  of  the  problem  are;  L  ~  5  cm,  c  ~ 
2.5  X  10^  cm/s,  u  ^  1  —  5  x  10®  s~^,  A  ~  0.05  —  0.25  cm.  Given  that 
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cell  sizes  must  be  on  the  order  of  A/4,  for  a  realistic  geometry  we  expect 
grid  dimensions  of  ~  100  —  500  pixels  in  each  dimension,  corresponding  to 
Np  ~  10"^  —  2.5  X  10^  in  2D  and  Np  ~  10®  —  1.25  x  10®  in  3D.  Furthermore, 
iterative  nonlinear  inversion  methods  may  require  many  evaluations  of  the 
forward  model  to  reach  convergence.  For  these  reasons,  development  of  highly 
efficient  methods  for  computation  of  the  forward  model  is  critical  to  the 
success  of  the  proposed  work.  A  significant  component  of  the  work  described 
in  the  following  report  has  centered  on  theoretical  analysis  of  the  scattering 
equation  and  computational  implementation  and  validation  of  the  forward 
model  in  both  2D  and  3D.  Our  model  is  capable  of  solving  the  full  scattering 
problem  on  a  desktop  computer  (2  GHz  PowerMac  G5)  on  a  1024  x  1024  grid 
in  2D  (see  Fig.  1)  or  a  100  x  100  x  100  grid  in  3D  in  approximately  30  minutes 
of  runtime,  and  has  been  structured  in  a  way  that  is  highly  amenable  to 
parallelization  in  anticipation  of  work  to  be  done  in  the  upcoming  months. 

2  Theoretical  Foundation 

2.1  The  Helmholtz  Problem 

We  begin  by  briefly  reviewing  the  basic  Helmholtz  problem: 

(A  +  Ko^(x))  ^p{x)  =  0  where  x  €  :  d  G  {2, 3},  (3) 

where  the  following  limit  holds  uniformly  in  all  directions: 

lim  (|x|^“V(x) -ffco)  =  0,  (4) 

|x|— >oo  \  / 

Here,  '0(x)  represents  the  total  pressure  field.  n(x)  is  the  refractive  index 
is  defined  by  n{x)  =  Co/c^(x),  where  cq  and  c(x)  are  sound  velocities  in  the 
ambient  medium  and  the  object,  respectively.  A  typical  scattering  geometry 
is  shown  in  Fig.  2. 

In  the  case  of  scattering  by  inhomogeneities,  it  is  customary  to  express 
the  Helmholtz  problem  in  terms  of  the  of  relative  change  in  the  refractive 
index.  Therefore,  we  define  a  quantity  7(x)  =  1  —  n(x)  =  1  —  ko/k(x),  where 
K  is  the  compressibility.  With  this  definition,  and  a  homogeneous  background 
medium,  we  observe  that  7(x)  has  compact  support  on  the  inhomogeneous 
scatterer  which  is  to  be  recovered.  In  a  non-attenuating  medium,  7(x)  is  a 
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pure  real  number,  while  an  attenuating  medium  has  a  complex  value:  7(x)  = 

7r(x)  +  i7i(x) 

The  formal  solution  of  the  Helmholtz  problem  may  be  expressed  as  a 
Lippmann- Schwinger  integral  equation: 

i/)(x)  =  7/;*"'^(x)  -kl  f  dx'G^^\ko\x  -  x'|)7(x')V'(x'),  (5) 

Jn 

where  fl  G  is  a  volume  bounding  the  scattering  object.  [1,  2]  The  incident 
scattered  wave  is  represented  by  a  plane  wave:  '^*"'^(x)  ==  where  the 

wavenumber  k^  =  |ko|  =  27r/A  and  A  is  the  wavelength.  'i/’(x)  represents 
the  total  field,  comprised  of  incident  and  scattered  components:  ^(x)  = 
^mc(x)  +.0«c(x),  G(0)(x|x,)  is  the  free  space  Green’s  function: 

G^°^(x|xs)  =  -^i/o(A:olx  -  x^l)  in  (6) 

and 

G(°)(x|x,)  =  inl^  (7) 

47r 

In  all  the  above  expressions,  x  is  the  field  or  observation  point  and  x^  is 
the  source  location.  Equation  5  is,  in  principle,  valid  for  all  x  G  R*^,  being 
applicable  both  within  and  outside  Q.  The  scattered  field,  ^®^(x),  at  the 
detector  position  X£)  can  thus  be  obtained: 

V/^(x£,)  =  -kl  [  dx'G(°)(A;o|x£,  -  x'|)7(x')'i/’(x')-  (8) 

Jn 

2.2  Discretization  of  the  Lippmann-Schwinger  Forward 
Problem 

Discretization  of  the  integral  equation  of  scattering  (Eq.  5)  was  performed  on 
a  regular  square  lattice,  though  preliminary  progress  toward  an  implemen¬ 
tation  on  a  hexagonal  lattice,  which  offers  the  benefit  of  increased  packing 
density  for  the  interpolation  of  the  Green’s  functions  (discussed  below),  has 
also  been  made.  Interpolation  of  the  Green’s  function  over  lattice-centered 
circular  elements  was  first  introduced  by  Richmond  for  the  2D  problem, [10] 
(see  Fig.  3)  and  has  been  extended  to  spherical  elements  in  3D  in  ths  work. 
For  circular  or  spherical  regions,  the  Green’s  function  can  be  analyticaly  in¬ 
tegrated  to  accurately  approximate  it  over  the  lattice,  significantly  speeding 
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up  the  initial  computation  of  the  discrete  Green’s  function  matrix,  the  accu¬ 
racy  of  our  integration  scheme  has  been  checked  against  full,  computationally 
expensive,  numerical  quadrature  of  G  over  the  individual  grid  cells.  In  addi¬ 
tion,  the  Richmond- type  discretization  makes  the  discrete  forward  problem 
Nystr6m-like,[ll]  which  allows  interpolating  the  nodal  values  of  the  solution 
to  arbitrary  points  in  the  lattice  thus  essentially  resulting  in  a  continuous 
solution. 

For  simplicity,  the  following  discussion  considers  the  2D  case  except  where 
there  is  a  significant  difference  in  3D.  Given  a  square  domain,  D,  of  side 
L  which  contains  the  scatterer,  and  Np  =  discrete  cells,  Eq.  5  can  be 
approximated  for  any  gridpoint,  {xjn,yn),  as: 

N  N  „ 

J/n)  “  ^  (^^rmyn)  ^0  /  dx'G^°)(A:o|x-x'|),  (9) 

p=l  q=l 


In  Eq.  (3),  a  =  L/N  is  the  lattice  spacing,  (xm,  yn)  =  a{m,  n)  is  the  coordi¬ 
nate  of  a  single  gridpoint,  {jtl>){xp^yq)  =  l{xp,yq)'4>{xp,yq)  and  the  integra¬ 
tion  domain,  F,  is  a  single  cell  {xmiVn)  l/n+i)-  Using  the  Richmond 

approach,  F  is  approximated  as  a  circle  of  area  equal  to  the  correspond¬ 
ing  cell  in  2D:  vq  =  :  M^,  and  as  a  sphere  of  volume  equal  to  the 

corresponding  cell  in  3D:  ro  =  (3/47r)^^^  : 

In  two  dimensions,  the  integral  in  Eq.  9  can  be  reduced  to 


dx'G^°^(kolx  —  x'l) 


^J,(koro)4^\hx)  x^x' 
^H[^\koro)  +  1  X  =  x'  ’ 


in  2D  and  to 


J  dx'G^°\ko\x-x'\) 


-i{koroyji{koro)h'i\kox)  x^x' 

1  —  —  ikoTo)  X  =  x'  ’ 


(10) 


(11) 


in  3D,  where  Jn{z),  H^\z),  jn{z),  and  h^\z)  are  the  standard  cylindrical 
and  spherical  Bessel  and  Hankel  functions.  [12] 


2.2.1  Evaluation  of  the  3D  Green’s  Function  Matrix 

The  task  here  is  to  evaluate  the  integral 

I=-kl  J  G(°)(A:o|x-x'|)da:'  (12) 
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of  the  free  space  Green’s  function  over  a  sphcric^al  resolution  element,  F'  = 
B{ro)  with  radius  ro,  which  separates  into  two  cases:  (1)  |x|  >  |x'|,  and  (2) 
|x|  =  |x'|  that  are  considered  separately  below. 


Off-Diagonal  Terms  (|x|  >  |x'|  :  x  ^  F'):  The  multipole  expansion  of 
(7(0)  gives 


oo  n 

G(°)(A:o|x-x'|)  =iA:o^  ^  a:„r„/i^^^{A:ox)P™(sin0)e“'”'^- 

71=0  7n=—n 


(13) 


where  x  =  (r,  9,  (p)  and  x'  =  (r',  6',  cp')  are  the  coordinates  of  the  field 
and  the  source  point,  respectively,  and  P^(x)  is  the  associated  Legendre 
polynomial. [12]  The  x'-integration  over  the  ball  yields 

f  jn{kox')Pn{sm  (^')e“*’"^'dx'  =  ^  (sin(A:oro)  -  (fco»'o)  cos{koro)) ,  (14) 

Jr'  ^0 

using  f  dzz’^jo{z)  =  sin(z)  —  zcos(2).  From  Eqs.  13  and  14,  the  off-diagonal 
elements  of  the  Green’s  function  matrix,  Qmh  are  found  to  be: 


=  -ih^oKkox)  (sin(A;oro)  -  (fcoro)  cos(fcoro)) 
=  -i{korofji{koro)h^o\kox). 


(15) 


Diagonal  Terms  (|x|  =  jx'j  :  x  G  F'):  In  order  to  calculate  the  diago¬ 
nal  elements,it  is  necessary  to  divide  a  resolution  element  into  two  separate 
regions  and  consider  the  integration  over  each  of  the  regions  separately.  The 
geometry  is  shown  in  Fig.  4.  In  region  I,  we  let  p  =  |x  —  £|,  where  i  is  the 
center  of  the  resolution  element  nearest  to  the  field  point  x.  The  calculation 
for  this  region  is  identical  to  that  for  the  off-diagonal  terms  computed  above, 
yielding: 

-klG^°'>  =  -ih^o\kop)  {sm{kop)  -  {kop)  cos{kop)) 

fi) 

=  -i{kopfji{kop)%  \kop). 

In  the  annular  region  II,  the  integral  is  similar,  noting  that  the  spherical 
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Hankel  function  may  be  written  as  h^\z)  —  jo{z)  +  iyo{z), 


pkoro 

=  -ijo{kop)  /  V) 

J  kop 

/  I'koro  pkoro  \ 

=  -ijo{kop)  f  /  drr‘^joir)+i  /  drr^yo(r)  J  • 

\J  kop  J  kop  J 

(17) 

This  equation 

is  evaluated  using  the  integrals 

J  dzz'^jo{z)  =  sin(z)  —  zcos(2), 

(18) 

and 

J  dzz^yoiz)  —  —  (cos(z)  -|-  2C0s(z)) , 

(19) 

giving 

(20) 

2.3  Solution  of  the  Lippmann-Schwinger  Equation  us¬ 
ing  the  CGFFT  Method 

Here  we  discuss  the  conjugate  gradient  fast  Fourier  transform  (CGFFT) 
method  as  applied  to  the  solution  of  Eq.  9.  If  we  attempt  to  solve  this  dis¬ 
cretized  Lippmann-Schwinger  equation  by  direct  matrix  inversion,  we  have: 

^  (21) 

where  A  is  a  Np  x  Np  matrix  having  the  values  of  7  along  the  diagonal: 
An  =  7i.  For  7  having  Np  =  elements,  the  operation  count  is 
making  direct  inversion  impractical  for  any  reasonable  sized  scatterer.  How¬ 
ever,  since  the  summation  in  Eq.  9  is  in  fact  a  d-dimensional  convolution 
it  may  be  evaluated  by  fast  Fourier  transform  (FFT)  in  0{N'^  log  N^)  op¬ 
erations,  leading  to  a  tremendous  speedup  for  non-trivial  grid  sizes.  [9]  This 
allows  us  to  formulate  the  forward  model  as: 

^  -  ^(0)  *  (^^),  (22) 

which  becomes,  after  Fourier  transformation, 

=  T  -  T  ,  (23) 
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so 


(24) 

We  solve  the  resulting  large  linear  system  iteratively  using  the  linear  con¬ 
jugate  gradient  approach.  [4,  5,  6,  7,  8]  The  Lippmann-Schwinger  operator 
is: 

Lip  = 'll:-  *  i'yijj).  (25) 

The  adjoint  Lippmann-Schwinger  operator,  can  be  calculated  from  the 
Hilbert  space  relation: 

{LiIj,(P)  =  {i/i,L^<P),  (26) 

where  the  inner  product,  {ip,  (p)  is 

{i/j,(p)  =  j  ip{x)cp^{x)dx,  (27) 

(p'^  denoting  the  complex  conjugate  of  (p.  Evaluating  the  inner  product  ex¬ 
pression  and  noting  that  G^'’^(x|x')  =  G^°^(x'|x)  yields  the  desired  adjoint 
operator,  namely 

(Z/'^0)(x)  =  0(x)  -  7(x)  f  G(°^^(A:o|x  -  x'|)0(x')dx'.  (28) 

Jn 

An  outline  of  our  CGFFT  routine  follows: 

•  Initial  step: 

Ro  =  LiPo  -  ij/^^ 

Pi  =  -L^Ro 

•  iterative  steps  (  k=l,2,..): 

(#"^LPfc)  ||L^Rfc_if 
""  lILPfcP  IlLPfcP 

ipk  =  ipk-i  +  afcPfc 
Rfc  =  L'i/'fc  —  =  Rfc_i  +  ojfcLPfc 

WL^Rkf 

||L-4p,_i||2 

Pfc+i  =  PkPk  — 
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The  residual  norm 


liRfcll  ||LV>fc-V/»-|| 

|j^inc|j  ||^*"C|| 

decreases  monotonically  as  the  algorithm  progresses,  and  is  an  useful  indica¬ 
tion  of  the  average  error  in  the  solution  after  k  steps. 

2.4  Analytical  Solutions  in  2D  and  3D 

In  order  to  verify  the  correctness  of  the  CGFFT  solution,  closed-form  bench¬ 
mark  solutions  were  computed  for  an  infinitely  long,  homogeneous  circular 
cylinder  of  radius  tq  in  2D,  and  for  a  homogeneous  sphere  of  radius  tq  in  3D. 

2.4.1  Scattering  From  an  Infinite,  Homogeneous  Circular  Cylin¬ 
der 

In  two  space  dimensions,  the  plane  wave  incident  field  can  be  rewritten  using 
the  Jacobi- Anger  expansion: 

OO 

V;*"^(r,0)=  ^  i"J„(A:r)cos(n0).  (29) 

n=-oo 

The  internal  field  is  then  given  by 

OO 

e)  =  Y^  hmJm{kor)  cos(m6'),  (30) 

m=0 


and  the  scattered  field  is: 


6>)  =  ^  amHl^\kor)  cos(m0).  (31) 

m=0 

We  let  ko  be  the  wavenumber  in  the  homogeneous  medium,  k  be  the  wavenum¬ 
ber  in  the  cylinder,  and  define  the  impedance  as  Z  —  k/kg.  The  coefficients 
am  and  bm  are  then  obtained  from  the  continuity  of  pressure  and  radial  ve¬ 
locities  across  the  surface  of  the  disk,  leading  to  the  following  expressions: 


—  fm  (^Jm{.kof^o)Jm{kr()^Z  Jm{kTQ^Jm 
bm  =  fm  (^Jm{koro)H^\koro)  -  H^\koro)Jm{koro)j 
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where 


frr 


Jmi^)  — 
Hm{z)  = 


j  1/ Aq  m  =  0 
\^2i^/Am  m^O 

Jm{kro)H^\koro)  -  H^\kQro)Jmikro)Z 
kim+liz')  Hrfi—l(^Z^. 


(33a) 

(33b) 

(33c) 

(33d) 


Here,  Hm\z)  =  Jmiz)  +  iYm(z)  is  the  first  kind  Hankel  function,  Jm{z)  is  the 
first  kind  Bessel  function,  and  Ym{z)  is  the  second  kind  Bessel  (or  Neumann) 
function. 


2.4.2  Scattering  From  a  Homogeneous  Sphere 

The  3D  case  closely  parallels  the  2D  case  discussed  above.  In  three  dimen¬ 
sions,  the  Jacobi-Anger  expansion  is: 

OO 

■ijr%r,  4>)  =  Y1  +  l)jn{k7  )Pn{sin  </>),  (34) 

n=0 

where  the  azimuthal  angle,  (p,  is  measured  with  respect  to  the  r^-plane.  The 
internal  field  is  then: 


OO 

7/;*'^*(r,  61)  =  J]]  bmjm{kor)Pm{sin  (p),  (35) 

m=0 


and  the  scattered  field  is: 


-0*^(r,  0)  =  a^j^)(A;or)F^(sin  (p).  (36) 

m=0 


Again,  the  coefficients  are  determined  by  requiring  continuity  of  pressure  and 
radial  velocity  at  the  interface,  giving: 

~  fm  (.?m(^O^o)jm(^^o)'^  jm(^kr())jiyi(^kQTo^'j 

(3T) 

bm  =  fm  (jm{koro)}f^\koro)  -  hk!i\koro)jm{koro)^ 
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where 


fr, 


1/ Ao  m  =  0 

2i'^{2m  +  1)/ Am  m  7^  0 


Am  =  jm{krQ)Um\hro)  -  h^mHkoro)jmikro)Z 


jmi^) 


((m  +  l)jm+l{z)  -  mjm-liz)) 


2m  +  1 

hmiz)  =  2m  +  I  -  m.hm-iiz)) . 


(38a) 

(38b) 

(38c) 

(38d) 


Here,  hm\z)  =  jm{z)  +  irimiz)  is  the  first  kind  spherical  Hankel  function, 
jm{z)  =  yJ-K j(^z)Jm-^\i‘i{.z')  is  the  first  kind  spherical  Bessel  function,  and 
Umiz)  =  \/Trl{2z)Ym+ll2{z)  is  the  second  kind  spherical  Bessel  (or  Neumann) 
function. 

2.5  Solution  of  the  Inverse  Scattering  Problem 

Our  ultimate  objective  is  the  efficient  and  accurate  reconstruction  of  the 
unknown  object  based  on  the  measured  scattered  field  (amplitude  and/or 
phase)  at  a  number  of  detector  locations  and  incident  frequencies.  Due  to 
the  incomplete  nature  of  the  measured  data  and  the  effect  of  measurement 
noise,  solution  of  such  inverse  problems  relies  on  the  use  of  ancillary  prior 
information,  typically  incorporated  as  an  ad  hoc  regularizing  term  in  $  which 
mediates  the  delicate  balance  between  the  agreement  between  the  simulated 
and  measured  scattered  fields,  ^  and  respectively,  and  the  ’’reasonable¬ 
ness”  of  the  reconstructed  object.  A  general  expression  for  the  objective 
function  is 

$(7)  =  $A*(7)  +  a$°(7),  (39) 

where  1Z  is  an  arbitrary  regularizing  functional  which  represents  the  concor¬ 
dance  between  the  reconstructed  object  and  our  prior  expectations  and  a  is 
a  regularizing  parameter  which  represents  the  weight  assigned  to  our  prior, 
relative  to  the  fit  quality,  Commonly,  is  chosen  to  be  the 
statistic  and  to  be  the  square  of  the  object’s  norm: 

^(7)  =  +  a||7|| 

=  ^(i>  -  -  r'") + o|j7p. 
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wIk^ic  E  is  the  covariance  matrix  of  the  measured  data.  For  discrete  obser¬ 
vations  with  uncorrelated  noise,  $(7)  becomes: 


^ob3 

*(7)  =  E 


m=l 


obs 


+  «  l7n 


n=l 


(41) 


where  am  is  the  standard  deviation  of  the  noise  for  the  m-th  detector.  In 
general,  the  objective  function  is  a  complicated  nonlinear  equation  which 
may  have  many  local  minima,  the  minimization  of  whiph  should  ideally  lead 
to  a  reconstruted  object  as  close  as  possible  to  the  true  object. 

As  pointed  out  in  the  introduction,  a  primary  concern  is  the  efficient 
computation  of  the  functional  (Prechet)  derivatives,  V-y4>,  of  the  objective 
function,  4>,  with  respect  to  the  object  to  be  reconstructed,  7.  Convention¬ 
ally,  this  gradient  would  be  evaluated  by  straightforward  differentiation: 


Np 

h  ^  ^  ^  h  7n*^m] 

n=l 


^obs 


(42) 


Obviously,  evaluation  of  the  full  Jacobian  matrix,  J,  at  every  iteration 
involves  significant  computational  expense  and  should  be  avoided.  This  can 
be  accomplished  by  the  computation  of  the  gradient  by  the  method  of  adjoint 
fields, [14,  15,  16]  which  is  briefly  described  here.  From  the  gradient  expression 
above,  the  total  variation  in  for  a  variation  J7  in  7  is  given  by: 


{^ohs 

m=l 


(43) 
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where  the  residual  Rm  is  defined  as  tjjrn  —  ^^nd  Si/j^  is  the  variation  in 
'ijjjn-  Now,  from  the  Lippmann-Schwinger  equation: 

=  (44) 

taken  over  the  computational  grid  Q.  Then 

Sipm^  G^^\xm\x)*{'Y5ip +  i;S^).  (45) 

Substituting  Eq.  45  for  S-ipm  into  Eq.  43  gives: 

(5$  =  fHe  +  ip5j) }  (46) 

in  Eq.  46  is  given  by 

^obs 

=  (47) 

m=l 

which  is  a  function  on  the  computational  grid.  Let  be  considered  as  the 
incident  field  in  the  Lipmann-Schwinger  equation.  Then: 

4>{x)  =  (x)  +  ,  70) ,  (48) 

where  0  is  the  total  field  corresponding  to  the  incident  field  0^*^^  Equation  48 
then  solves  0  as: 

0=  [/-^(o)A.^]-yo).  (49) 

But  from  the  Lipmann-Schwinger  equation,  we  also  have: 

50  =  [/  -  g^A.,]- V,  (50) 

which  leads  to,  after  some  algebra, 

S(p  =  Re  I  dx{(f){x)'iJj(x))S'y{x) 

Ju 

—  Re  /  dx{V'y^)S'y{x), 

Ja 

from  which 

V.y$(x)  =  i?e(0(x)0(x)),  X  E  ft  (51) 


13 


3  Numerical  Results 

Here  we  compare  the  results  of  our  CGFFT  solutions  of  the  Lippmann- 
Schwinger  equation  for  both  and  with  the  results  of  Born 

approximation  calculations  in  2D  and  with  the  corresponding  exact  analyt¬ 
ical  expressions  derived  in  the  previous  section  in  both  2D  and  3D.  We  also 
demonstrate  simulations  for  scattering  from  a  square  object  in  2D  for  vari¬ 
ous  incidence  angles.  Finally,  we  present  preliminary  results  of  performance 
profiling  and  discuss  the  significance  for  the  inverse  problem. 

3.1  Comparison  with  Born  Approximation  and  Exact 
Results 

Figure  5  shows  results  for  the  absolute  magnitude  of  the  scattered  field,  along 
with  the  real  and  imaginary  components,  from  an  infinite  circular  cylinder 
in  2D  for  two  different  parameter  regimes.  Born  approximation  results  are 
shown  by  the  blue  dots,  the  exact  calculation  by  red  dots,  and  the  CGFFT 
computation  by  black  circles.  Plots  on  the  left  side  of  the  figure  are  for  a 
parameter  regime  well  within  the  realm  of  applicability  for  the  Born  approxi¬ 
mation,  as  is  evidenced  by  the  extremely  high  degree  of  concordance  between 
all  three  computations.  Where  any  discrepancy  at  all  is  discernible,  it  is  the 
Born  data  which  deviate  slightly  from  the  exact  and  CGFFT  results  which 
are  essentially  indistinguishable.  The  right  half  of  the  figure  presents  a  much 
more  challenging  test  case  in  a  parameter  regime  where  the  Born  approxi¬ 
mation  is  clearly  invalid.  Here,  the  plots  confirm  our  expectations;  the  Born 
data  deviate  dramatically  from  the  exact  solution.  However,  as  before,  the 
results  of  the  CGFFT  solution  are  virtually  identical  to  the  exact  values  for 
the  scattered  field.  In  Fig.  6,  we  show  similar  results  for  scattering  into  a 
ring  of  detectors  oriented  perpendicular  to  the  r0-plane  by  a  spherical  ob¬ 
jects.  As  before,  the  exact  and  CGFFT  computations  agree  perfectly  in  the 
Born  limit.  In  the  non-Born  regime,  we  begin  to  observe  noticeable  devia¬ 
tions  between  the  exact  and  CGFFT  solutions;  however,  as  we  have  verified 
with  numerical  experimentation,  these  deviations  stem  from  the  necessarily 
much  coarser  discretization  of  the  object  in  3D.  Despite  these  small  discrep¬ 
ancies,  it  is  clear  that  the  CGFFT  method  is  capturing  the  essential  physics 
of  scattering  in  the  3D  case  as  well  as  in  2D. 

To  more  accurately  quantify  the  level  of  accuracy  we  can  expect  from 
our  CGFFT  algorithm,  we  present  a  comparison  of  between  CGFFT 
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and  the  exact  solution  for  the  infinite  cylinder  and  sphere  in  Fig.  6.  For  the 
circle,  computed  on  a  512  x  512  lattice,  we  find  essentially  perfect  agreement 
between  the  two,  with  a  maximum  discrepancy  of  less  than  0. 1  percent  and  an 
average  deviation  of  less  than  0.02  percent.  Again,  due  to  the  more  prominent 
discretization  effects  (both  the  circular  and  spherical  objects  were  smoothed 
with  a  0.5  pixel  Gaussian  blur  filter  to  reduce  scattering  artifacts  arising  from 
pixelation  at  the  object  boundary)  arising  from  coarser  grid  size  in  the  3D 
case,  the  error  is  enhanced  relative  to  the  circular  object.  Despite  that  fact, 
the  absolute  value  is  bounded  at  less  than  1  percent,  and  averages  less  than 
0.4  percent.  In  addition,  the  exact  and  CGFFT  internal  fields  are  clearly 
qualitatively  quite  similar.  We  demonstrate  the  ability  to  simulate  more 
complex  objects  in  Fig.  8,  which  demonstrates  changes  in  both  the  internal 
and  scattered  fields  with  rotation  of  the  angle  of  the  incident  radiation,  as 
would  be  expected  for  an  anisotropic  object. 

3.2  Performance  Profiling 

In  order  to  assess  the  performance  of  our  forward  model  and  identify  the 
presence  of  performance  bottlenecks,  we  have  investigated  timing  behavior 
in  a  number  of  cases.  In  Fig.  9,  we  present  the  dependence  of  execution 
time  for  fixed  7  and  A  at  various  lattice  sizes.  Plotting  total  execution  time 
against  Np,  the  number  of  gridpoints,  reveals  behavior  consistent  with  the 
Np  log  Np  expected  for  an  algorithm  which  is  dominated  by  the  cost  of  per¬ 
forming  forward  and  inverse  FFTs,  with  a  few  notable  discrepancies.  The 
clear  appearance  of  outlier  data  points  with  substantially  elevated  execution 
time  is  connected  with  FFTs  on  data  sets  which  have  a  large  prime  factor; 
this  is  a  well  known  effect  in  FFT  optimization  and  emphasizes  both  that  our 
algorithm  is  FFT-performance  limited  and  that  it  is  of  great  importance  to 
select  grid  sizes  with  small  prime  factors.  Another  interesting  outlier  is  the 
significant  performance  penalty  associated  with  the  128  x  128  FFT,  which 
is  50%  slower  than  would  be  expected  from  the  performance  curve.  In  this 
case,  it  appears  that  memory  access  collisions  in  the  fast  cache  memory  are 
causing  an  anomalous  slowdown. 

Two  additional  effects  are  observed  to  have  a  significant  impact  on  per¬ 
formance,  both  relating  to  the  convergence  rate  of  the  conjugate  gradient 
algorithm  and  understandable  as  arising  from  progressively  greater  devia¬ 
tion  from  the  Born  approximation  which  is  used  as  a  starting  point  in  our 
conjugate  gradient  algorithm.  The  first  is  seen  with  increasing  7,  which  cor- 
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responds  to  progressively  stronger  scattering  from  inhomogeneities  with  re¬ 
fractive  properties  that  differ  more  from  the  background  medium.  Figure  10 
shows  a  quadratic  dependence  on  computation  time  with  increasing  values 
of  7.  Similarly,  in  the  second  case  of  decreasing  wavelength,  A,  we  expect 
to  see  increases  in  computation  time  as  the  spatial  scale  of  structure  in  the 
internal  field  increases.  Figure  11  plots  the  A-dependence  for  three  different 
2D  lattice  sizes,  showing  consistent  power  law  behavior,  with  execution  time 
varying  approximately  as  1/A^-^. 

4  Key  Research  Accomplishments 

•  Complete  theoretical  analysis  of  problem  formulation. 

•  Rapid  adjoint  computation  of  gradient  terms  validated  against  finite 
difference  results. 

•  Implementation  of  efficient  CGFFT  algorithm  for  2D  and  3D  forward 
modeling. 

•  Development  of  approximate  and  exact  analytical  solutions  in  2D  and 
3D  and  validation  of  CGFFT  algorithm  against  them. 

•  Approximate  expression  for  Green’s  function  validated  against  result 
of  full  numerical  quadrature. 

•  Implementation  of  nonlinear  conjugate  gradient  method  for  solution  of 
inverse  problem  in  2D  and  3D  complete. 

•  Initial  investigation  of  various  forms  of  regularizing  functional  and  var¬ 
ious  approaches  to  optimizing  the  convergence  of  minimization  of  the 
objective  function. 

5  Reportable  Outcomes 

•  Abstract  and  oral  presentation  at  the  16th  Annual  UCAIR  Symposium, 
Park  City  UT,  October  2004. 
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•  Abstract  submitted  to  2005  Fully  Three-Dimensional  Image  Recon¬ 
struction  Meeting  in  Radiology  and  Nuclear  Medicine,  Salt  Lake  City 
UT,  July  2005. 

•  Abstract  submitted  to  Applied  Inverse  Problems  2005,  Cirencester  UK, 
July  2005. 

6  Conclusions  and  Future  Work 

From  the  discussion  above,  it  is  apparent  that  a  sufficiently  accurate  rep¬ 
resentation  of  the  prostate  will  require  a  resolution  for  which  the  compu¬ 
tational  burden  will  exceed  reasonable  limits  for  a  single  processor  system. 
Fortunately,  much  of  the  CGFFT  computation  is  readily  parallelized  to  take 
advantage  of  multiprocessor  hardware  and/or  the  high-performance  comput¬ 
ing  cluster  systems  that  are  rapidly  growing  in  popularity  and  availability. 
Nevertheless,  while  such  approaches  are  likely  to  be  relevant  to  this  work, 
it  remains  critical  to  investigate  all  avenues  to  improving  the  intrinsic  algo¬ 
rithm  performance  before  resorting  to  brute  force.  For  the  upcoming  year, 
we  will  work  in  a  number  of  areas: 


6.1  Improving  Forward  Model  Performance 

•  preconditioning  strategies  for  acceleration  of  conjugate  gradient  con¬ 
vergence 

•  implementation  of  the  Non-Uniform  FFT  to  enable  the  use  of  unevenly 
spaced  grids 

•  predictor-corrector  methods  for  minimizing  evaluations  of  full  forward 
model 

•  parallelization  of  forward  model  algorithm  and  testing  on  cluster  system 

•  extension  of  theory  to  accommodate  point  source  incident  waves  and 
transurethral  geometry 

•  extension  of  theory  and  model  to  incorporate  mass  density  terms  and 
attenuation 
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6.2  The  Inverse  Problem 


We  have  implemented  and  performed  preliminary  tests  on  the  adjoint  gra¬ 
dient  inversion  algorithm  in  both  2D  and  3D,  demonstrating  the  ability  to 
reconstruct  simple  geometries  in  fairly  low-resolution  cases.  However,  there 
is  much  improvement  to  be  made  in  the  inversion  algorithm.  Among  the 
strategies  which  we  intend  to  employ  over  the  next  year  are: 

•  investigation  of  alternative  and/or  hybrid  minimization  strategies  such 
as  simulated  annealing  and  genetic  algorithms  to  be  used  in  conjunction 
with  nonlinear  conjugate  gradient  methods  to  maximize  convergence 
rate. 

•  study  of  different  regularizing  functionals  including  total  variation, 
maximum  entropy,  laplacian,  and  others  and  investigation  of  their  im¬ 
pact  on  inverse  problem  convergence  properties. 

•  analytical  and  numerical  solution  of  the  variational  problem  in  3D. 

•  parallelization  of  inversion  algorithm 

•  development  of  a  realistic  prostate  phantom  using  MRI  and  conven¬ 
tional  ultrasound  data  in  conjunction  with  anatomic  information 
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Figures 
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Figure  1:  Internal  field  amplitude  eomputed  for  a  1024  x  1024  circular  scat- 
terer  object  with  7  =  0.1  and  \/L  —  1/28  for  a  single  plane  wave  incident 
along  the  a;-axis. 


Homogeneous  background  :c{x)  = 


Figure  2:  Schematic  illustration  of  typical  geometry  of  an  inhomogeneity 
scattering  problem.  An  incident  plane  wave  traveling  along  the  x-axis  im¬ 
pinges  on  an  object  embedded  in  the  homogeneous  background  medium, 
which  is  bounded  by  the  problem  domain,  fl.  Outgoing  scattered  waves 
travel  in  all  directions  from  the  object,  and  are  measured  at  detector  posi¬ 
tions  arrayed  outside  Q. 
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Figure  3:  Schematic  illustration  of  the  discrete  square  lattice  in  2D.  The 
circle  (extended  to  a  sphere  in  3D)  of  radius  ro,  represents  the  region  over 
which  7  is  considered  homogeneous  and  the  Green’s  function  is  analytically 
integrated.  The  radius  shown  is  smaller  than  the  actual  radius  for  clarity; 
the  actual  value  of  ro  corresponds  to  an  area-preserving  integration.  L  is  the 
physical  size  of  the  domain,  Q,  with  a  lattice  spacing  of  a  =  L/N.  Element 
centers  are  nodes  for  the  basis  functions  in  a  Galerkin-type  projection. 
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Figure  5:  Comparison  of  ipse  for  CGFFT  (black  circles),  exact  solution  (red 
dots) ,  and  the  Born  approximation  (blue  dots)  for  scattering  from  an  infinite 
cylinder,  represented  on  a  128  x  128  grid.  The  upper  panels  plot  \ipsc\,  the 
middle  plot  fHe{'0sc},  and  the  bottom  panels  plot  3m{ijjsc}-  Panels  on  the 
left  are  for  L  =  li  cm,  A  =  28  cm,  and  7  =  0.01,  while  panels  on  the  right 
are  for  L  —  M  cm,  A  =  3.5  cm,  and  7  =  0.3. 


Figure  6:  Comparison  of  ijjsc  for  CGFFT  (black  circles),  exact  solution  (red 
dots),  and  the  Born  approximation  (blue  dots)  for  scattering  from  a  sphere, 
represented  on  a  32  x  32  x  32  grid.  The  upper  panels  plot  \ijjsc\,  the  middle 
plot  and  the  bottom  panels  plot  Panels  on  the  left  are 

for  L  =  14  cm,  A  =  28  cm,  and  7  =  0.01,  while  panels  on  the  right  are  for 
L  =  14  cm,  A  =  3.5  cm,  and  7  =  0.3. 
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Figure  7:  Comparison  of  ijjint  for  CGFFT  and  the  exact  solution  for  scattering 
from  an  infinite)  cylinder,  represented  on  a  512  x  512  grid  in  the  left  panels 
and  from  a  sphere,  represented  on  a  a  64  x  64  x  64  grid  in  the  right  panels. 
The  upper  panels  plot  \ipint\  from  the  CGFFT  calculations,  the  middle  plot 
I'lpintl  for  the  exact  solution,  and  the  bottom  panels  plot  the  relative  error, 
{\'ipint,CG\  -  \'ipint,ex\) /\i’int,ex\-  In  both  cases,  the  forward  model  was  run  with 
L  —  14  cm,  A  =  7  cm,  and  7  =  0.3. 


Figure  8:  Internal  (upper  panels)  and  scattered  field  (lower  panels)  ampli¬ 
tudes  computed  for  a  128  x  128  square  object  with  7  =  0.3  and  A/L  =  1/2 
for  three  plane  waves  having  incidence  angles  of  0,  7r/8,  and  7r/4. 


2.5 


Forward  model  performance  :  dependence 


Figure  9:  Scaling  of  forward  model  execution  time  with  number  of  gridpoints 
in  the  scattering  object.  While  some  expected  deviation  is  seen  for  smaller 
sizes,  the  overall  trend  can  be  seen  to  be  dominated  by  the  Np  log  Np  contri¬ 
bution  from  the  FFT.  Furthermore,  the  outlier  points  occur  for  FFT  lengths 
whose  greatest  prime  factor  is  .F  G  {13,17,19,23,29,31},  leading  to  well- 
known  inefficiencies  in  FFT  implementations,  and  for  an  FFT  length  of  128, 
which  presumably  causes  cache  collisions  in  the  2D  case  presented  here. 


Figure  10:  Scaling  of  normalized  forward  model  execution  time,  i„orm  = 
tex/NplogNp,  with  7,  demonstrating  a  quadratic  dependence  of  algorithm 
convergence  with  increasing  deviation  from  Born-regime  scattering. 


Forward  model  performance :  X  dependence 
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Figure  11:  Scaling  of  normalized  forward  model  execution  time,  tnorm  = 
tex/NplogNp,  with  A  for  three  values  of  N  :  {64,80,96}.  Execution  time 
is  shown  to  obey  a  power  law  behavior,  independent  of  Np.  Regression  to 
the  linear  portion  of  the  curves  leads  to  a  power-law  relation  of  the  form 

inorm  ^ 


